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PROGRESSIVE WAVE EXPANSIONS AND OPEN 
BOUNDARY PROBLEMS 

T HAGSTROM* AND S I HARIHARAN* 

Abstract. In this paper we construct progressive wave expansions and asymptotic 
boundary conditions for wave-like equations in extenor domains, including applications 
to electromagnetics, compressible flows and aero- acoustics The development of the con- 
ditions will be discussed in two parts The first part will include denvations of asymptotic 
conditions based on the well-known progressive wave expansions for the two-dimensional 
wave equations A key feature m the denvations is that the resulting family of bound- 
ary conditions involve a single denvative in the direction normal to the open boundary 
These conditions are easy to implement and an application in electromagnetics will be 
presented The second part of the paper will discuss the theory for hyperbolic systems 
in two dimensions Here, the focus will be to obtain the expansions in a general way and 
to use them to derive a class of boundary conditions that mvolve only time derivatives 
or time and tangential derivatives Maxwell* s equations and the compressible Euler 
equations are used as examples Simulations with the linearized Euler equations sure 
presented to validate the theory 

Key words. Progressive wave expansions, boundary conditions, Maxwell’s Equa- 
tions, Euler Equations, Numerical Simulations 

AMS(MOS) subject classifications. 65M99, 35B40 


1. Introduction* Exterior problems are commonly posed for wave- 
like equations, and their numerical solution leads to the problem of open 
boundary conditions. We discuss both isotropic and nonisotropic cases as 
they arise in electromagnetics and fluid dynamics. These equations in- 
clude first order hyperbolic systems such as Maxwell’s equations, the Euler 
equations of compressible flows, or the linearized Euler equations, as well 
as second order reduced forms as appropriate Many work studies of this 
problem have appeared m the recent literature and we won’t try to list them 
all There are fundamentally two different, though related, approaches that 
have usually been taken One is the use of high frequency asymptotics such 
as the geometrical optics approximation. The other is based on the far field 
structure of the solution. (For a third approach based on the direct ap- 
proximation of the exact condition, see [5].) Progressive wave expansions 
were used as a tool to construct far field boundary conditions as early as 
the time of Sommerfeld In the modern computational pomt of view, they 
were put in use for the first time by Kriegsmann and Morawetz [8] Since 
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then there have been many variations to this approach. For example, ex- 
tensions to anisotropic propagation were first attempted by Bayliss and 
Turkel [2], and a generalization to the case of anisotropic wave equations in 
two and three dimensions was proved by the authors [7] Issues include the 
construction of the expansion for general systems, their use to construct 
stable boundary conditions of minimum order, and, finally, their practi- 
cal implementation. Higher order conditions no matter which approach is 
used, are typically more complicated than the partial differential equation 
one starts with In particular, they tend to have higher order derivatives m 
the direction of the propagation. To avoid this problem, often the partial 
differential equation itself is used. Such a procedure is not known in general 
for problems governed by first order hyperbolic systems. Here we provide a 
systematic way of dealing with this issue using progressive wave expansions. 
Our attention focuses on first order systems, namely, Maxwell’s equations 
and the linearized Euler equations To motivate the central ideas, we first 
consider the second order wave equation with the emphasis on progressive 
wave solutions. 

2. Second order wave equation. As mentioned above, the goal 
here is to treat the problem of boundary conditions without having higher 
order normal derivatives. To illustrate the underlying procedure, let us 
consider the problem governed by the wave equation m two dimensions. 
We wish to construct the progressive wave solutions to this equation and 
exploit their structure to prescribe asymptotic boundary conditions. The 
equation written in cylindrical coordinates takes the form 

(2.1) Utt = Urr + -Ur + 

t r* 

We look for solutions that are periodic in the angular direction as follows 

oo 

(2.2) u(r,9,t) = ^v n (r, t) a„ (0) 

n= 0 


where On(0) = A n cos nO + B n sin nO. Substituting (2.2) in (2.1), we obtain 

1 n 2 

(2.3) Vn y tt = V ny rr + ~V n ,r - ^jV n 


Following Friedlander [4], we construct solutions of (2.3) in the form 

(2.4) ”«(r,<) = f 


J=0 


Substitution of (2.4) in (2 3) results m the following recurrence relations 

f W , t X _ (j+?) 2 -n 2 

f 3+ i(t r)-- 2(j + 1) 


(2.5) 
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The goal is to examine the effect of the recurrence relations on constructions 
of asymptotic boundary conditions. First, we observe that substitution of 
(2.4) in (2.2) yields the following formal representation of the solution 

00 00 f n (t-r) 

(2.6) w(r, 0, t) = 2 X 

n=0 j=0 T 2 


Mampulations of this senes, particularly increasing the order of the decay 
rate for boundary conditions, have been proposed by many authors (e.g. 
Bayliss and Turkel [1]). In fact, a different form of (2 6) has been used for 
these constructions, which will not be discussed here. We define a “basic 
boundary operator” from (2.6) as follows: 

„ d d 1 
B = di + fr + Tr 

It is immediately venfied from (2 6) that 



( 2 . 8 ) 


n=0 j=l 


r j+§ 


Direct approximation of (2 8) is the radiation condition Bu = 0, a pop- 
ular condition in the literature noted by many researchers (e g. Bayliss 
and Turkel [1], Engquist and Majda [3]). Asymptotic accuracy of such 
a condition is 0(r -5 / 2 ), which is evident from (2.8). Higher order condi- 
tions in general require higher order normal derivatives or derivatives in 
the dominant direction of propagation This may not be a desirable fea- 
ture numerically, particularly for nonlinear generalizations. Here we obtain 
higher order conditions that involve Bu, u, u$$, and their time derivatives 
on the artificial boundary. We begin with the construction of higher order 
conditions by differentiating (2 8). This yields: 


(2.9) 


4#“= -X a „(<?)X 


dt 


n=0 


J=U 


r J+f 


Noting that the inner summation may be written in the form 


oo 


E 


0 + !)/”+! (* - r ) 
r j+§ 


and using the recurrence relation (2 5) yields: 




n=0 
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A simple manipulation of the right hand side yields: 

oo 


( 2 . 10 ) 


d 1 1 6Pu _ 1 V" (ns. + 

BU &r* U 2r*d9 2 ~ T J+f f > 

n= 0 j=l 


The highlight here is the observation 

00 °° .ff(t-r) 


( 2 . 11 ) 


«** = -51 »*«••(*) 2 ; 


n=0 


3=0 


W + 2 


Note that the asymptotic accuracy of the candidate boundary condition 
(2 10) is increased further to 0 (r” 2 ) Let 


( 2 . 12 ) 

Then (2 10) becomes 
(2.13) 


d u u ee 

** “ a*— «? " -? 




n =0 


3=0 


This form again suggests the use of the recurrence relations (2 5) by differ- 
entiating the equation with respect to time Doing so, we obtain 


(2 14) 




n =0 


3=0 


which is equivalent to 


<*■»&*• + - - 1 E ■ - <*> E 'fr+ff-’* 1 /; 


n =0 


J =1 


We note that — n 2 translates into the second tangential denvative 
Defining 


(2.16) 


U+ £* U+ ^ 3 U " 


it is clear that equation (2.15) yields a one asymptotic order higher bound- 
ary condition (to 0(r“f)). Moreover, noting 


(2.17) 


b = -1 (j + - 1 - 0 - + §) ’ - ”’ 1 " 

n =0 j =0 




J + l 


and applying the time derivative again, the process becomes clear and it 
yields 

„ d 25 13 1 1 x 

(2 18) S 3 « .= jB 2 u - - —u ee - —ueeee = 0(-^) 
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Remark 1: As far as numerical implementation of these conditions are 
concerned, one may consider a sequence of operations to update u at the 
current time: 


(2 19) 


1 

U t + U r + T-tl = Z 


( 2 . 20 ) 

(2.21) 


z *-fr2 u+ ^ u »* =v 


t ’ t+ 8^ U+ 2^ UM = W 


25 13 1 

(2 ' 22) " 128 ^“ - T^ U " ~ & Utm = 0 

The above sequence of equations (which provides a boundary condition 
asymptotically accurate to 0(r“^ )) as a system of first order equations to 
march in time. (At the time this article was written one of the students 
of the second author has implemented such a procedure and obtained the 
indicated asymptotic improvement. The details will appear elsewhere). 
Remark 2: The procedure above coincides with the high frequency ap- 
proximations of the exact condition in the radially symmetric case In the 
Laplace transform domain, the exact operator has the form (see [6]) 


(2.23) 


B c u = 


1 stKq(st) 
r K 0 (st p 


Where Kq ( z ) is the modified Bessel function of order 0. Moreover we find 
as sr — ► oo 


(2 24) = + 37 


+ 0((«r) )) 


The Laplace transform of the derived operators coincides with the large sr 
approximations of the exact boundary operator B e We can, then, interpret 
the expansions both as a long-range and as a high-frequency approximation. 

We also note that the Fourier transform of the operator B\u coincides 
with the second order operator proposed by Knegsmann et al. [9] in con- 
junction with on surface radiation conditions. As an example we consider 
the computation of the surface current calculation in electromagnetic scat- 
tering. Let T be the boundary of a perfect conductor Then the magnitude 
of the total current is given by the formula (see [9]) 
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where u 3 is the scattered field, and u xnc is the known incident field For 
perfect conductors u 3 = — u tnc on the boundary of the scatterer T The 
principle of the on surface boundary procedure consists of bringing the far 
field boundary exactly on the interface of the scatterer The advantage is 
rather clear. Since the total current is a functional of the normal derivative 
of the scattered field and the radiation boundary operators on the surface 
directly express the normal derivatives m terms of the mcident field We 
note that m the formula for the surface current k is the wave number which 
arose from the Fourier transform of the wave equation We list the Fourier 
transform of the operators denved m our theory They are: 

Condition 1 

(2 26) 

Condition 2. 

(2 27) - tk(-tku + u r + ^ u ) = ^ + 

Condition 3. 

(2.2$Q -*&) 2 (—tku + u r + ^ u ) = -*fc(g ^2 w + - ^3 u »» 

The first two operators are used m [9], and the third one is, so far as we 
know, new A plane wave mcident upon a unit cylinder is considered for 
the calculation of J and results are shown in Figures 2.1 ( k = 5) and 2 2 
( k = 2) respectively. The incident field has the specific form u tnc = e xkrcosB 

Remark 3.: For anisotropic equations, such as convective wave equation, 
an analogous procedure may be denved. The use of the resulting conditions 
are more pertinent to systems of equations such as the Unearned Euler 
equation This is discussed m section 4. 


- tku + u r 4- —u = 0, 
2r 



Surface Current Surface Current 
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Surface Current J (k=5 0) 


Comparisons with ’Exact SofcAon’ 



FlG 2 1 Comparison of results with exact solution , k = 5 


Surface Current J (k=2 0) 


Companson with Exact Solution 



FlG 2 2 Comparison of results with, exact solution, k = 2 
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3. First order hyperbolic systems - Isotropic case. Here our 
focus is to extend the ideas to systems of first order equations The pro- 
gressive wave expansions may be carried out directly in the time domain 
as we did for the second order wave equations or m the Laplace trans- 
form domain. In this section we present the construction using the Laplace 
transform. The direct approach is illustrated in Section 4. 

Maxwell’s equations offer an interesting example of an isotropic sys- 
tem Here we confine our attention to Ttansverse Magnetic (TM) fields for 
simplicity. The full field equations are: 

(3 1) div eE = dtv /zH = 0, 

(3.2) curl H = 


(3.3) curl E = 

We shall consider TM fields as follows: 

(3 4) E = E(x,y,t) k 


dE 

'dt ’ 

du 


-M 


at 


(3 5) H = fTi(*,y,t)i + H 2 (x,y,t) j. 


Equation (3.4) indicates that the electric field propagates m the direction 
perpendicular to the x — y plane and is transverse to the magnetic field 
Under these assumptions, equations (3.2) and (3.3) become 


(3 6) 

and 

(3.7) 


dH 2 _ dHi _ dE_ 
dx dy € dt 


dE_ 

dy 


. dE. (dH t . ^ dH 2 A 


respectively Rearranging equations (3 6) and (3 7), we obtain the following 
system 


dE _ 

1 

dH 2 

dE{ 

dt ~ 

€ 

dx 

dy . 


dHi _ 1 dE 

dt dy 


dH2 _ 1 dE 
dt (i dx 
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This can be put in the conventional form: 


where 


u t = Au,+Bu y 


0 0 1 
0 0 0 
I o 0 


0 0 
-i 0 0 

0 0 0 


and where u = (E, Hi,H 2 ) t Converting to polar coordinates we obtain. 


« u r 


-ism0 

- COS0 
A* 


— ^sin0 ^cos0 \ 

0 0 I _ 

0 0 / * 


+ - I -^cos 6 

T V 


-icos0 — 7 sin# \ ( E 

S 2 Hu 


We take the Laplace transform of (3.9). With the change of variable f = rs 
we have 


(3.10) 


where 


= (*& + ?©£) A 
V b 2 


— TSin 0 t cos# 


R - I -isrntf 

V Jl cose 


©= j ^ cos 0 

V -^sind 


— t cos 6 — 7 sin 0 


We seek an expansion of solutions of (3 10) in the form 


ft = V~( a+ " b+ - ) 

#2 / 
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We note that this form is similar to Fnedlander’s form that applies to 
the second order wave equation in the Laplace transform domain Also, we 
have introduced a decay rate constant a which turns out in two dimensions 
to equal as expected Substituting (3.11) into (3.10) yields to leading 
order: 

(3 12) A = I + gR + g f ® 

and 

(3 13) Aa=0 

For this requirement to be true, clearly it must follow that* 

C l — ^ (g sin 0 + g f cos 0) 

-i( 3 sin 0 + g'cos6) 1 

±(gcos0 - g'sinff) 0 


cos 6 - 
0 
1 


= (1 — — ((psin# + < 7 'cos 0 ) 2 + ( gcosO — g' sw0) 2 )) 

€fl 


(3.14) 

The roots of equation (3.14) are: 


= 1 -j-y+w?) 


g = ± v/i/I, y/ejl cos(0 + <£), 

<f> arbitrary. For waves propagating to infinity in all directions we choose 
g = y/e^L as the allowoable root With this value of p, the matrix A becomes 


A = 


( 1 — y^sinfl ^cos6 ^ 

1 0 

Jj-cosO 0 1 


whose right nullvector is 


a — 


( i 

fisinO 

-Jocose 


ai(0) 


with left nullvector 


= (1, ^j-sinfl, — ^~^cos6) 


g r sin 6) 
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The next order terms m the asymptotic expansion yield the following rela- 
tion: 


(3 15 ) 

and a is determined by 
( 3 . 16 ) 


Ab = -aRa + e^, 


a = 


i T eff 

l T Ra ’ 


ie., by requiring (1, Ab) = 0 

Noting the following calculations: 


Ra = a\(0) 



da 

do 


( 

dai 

lo 

\ 


1 

\ft sm9 

-yj 'Jcosff 


\ 


0 

\ 

4 - &i 

\ 

/jc os 6 


/ 

l \ 

f^sine 

) 


da da i 

& d$ ~ He 



l T Ra = a 1 (0)(-J ? ), l T e| = oaW(^=), 


W 


W 


it follows that 


Choose 



and use the last two equations to obtain 
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, . 1 .. 1 a dai 

h = (- — cos^jai smO—. 

2 n (i ov 

Substituting a and b in (3.11) we obtain 

_ e _f v^ BE _ e" f v^ 3ai 
E ~ fi ai ’ dO ~ fi dO ’ 


fr /T . - 1 .sin0 - .di? , 

I7i = ,-smOE -\ (— — E - cos 6-^-), 

V M hts k 2 dO’ 

fr /T - 1 . — cos# - n dE . 

H 2 = -cosOE+ ( — - — .E-sm0— ). 

Y /i 2 <90 


Multiplying through by s and taking the inverse transform finally yields 


(3 17) 
(3 18) 


dH\ , IT .d sin 6 

-ar = ( VJ“ # a + 2?r 


dH 2 

dt 


.IT n d cos 0 

= -<V^ * + ^ 


cos# d 
\ir dO 


)E , 


sin 0 9 . _ 

~jtTde^ E 


As only one boundary condition is required, we convert these into a 
single condition Multiplying (3.17) sm0 and subtracting (3.18) multiplied 
by cos0, we obtain our final form* 


(3.19) 


8H 2 

— — sm 0 — cos 0 - 

dt dt 




This construction is easily extended to higher order, though we have 
not devised a unified approach to the implementation of the higher order 
conditons. 


4. The linearized Euler equations - An anisotropic example. 
The construction of asymptotic boundary conditions for the linearized and 
nonlinear compressible Euler equations is also of interest, particularly for 
apphcations m aeroacoustics. In this section, we construct the expansions 
in the time domain directly. Again, the system takes the form: 

(4.1) u t = A n x + B u y 

where A and B are constant matrices In cylindrical coordinates we have 

U t = R Ur + -T u$, 
r 


(4.2) 
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where 


(4.3) 

R = Acos0 + Bsin0, 

(4.4) 

T = —A sm^ + B cos 0 

A far field asymptotic expansion may be constructed in the following form 
for the solution vector u: 

(4.5) 

u ^Ut-rg(0)) 

U - 2^ _„+i ®»» 

n=0 

where the scalar function g(9) and the vectors an(0) axe to be determined 
The function /o is analogous to the radiation function discussed m [4] The 
other functions are recursively determined by substitution of the expansion 
into equation (4.2). The 0(£) terms yield: 

(4.6) 

Cao — 0, 

where C 
det( C) = 
yields 

= I + </(0)R + g f (0) T For ao to be nontrivial one must have 
0, yielding an ’eikonal function’ g(0). The next order correction 

(4.7) 

/oC ai = -/(iRao + T^ 


This imposes a necessary restriction that fo = f{ In general, it follows 
/n-i = /A? 71 > 1* At this point, we turn to the lsentropic, linearized, 
compressible, Euler equations to illustrate the actual calculations involved 
in solving these algebraic problems For a uniform base flow in the x 
direction they are* 


(4.8) 

, d . _ d . du dv 

{ di +M di )p+ di + d^ = Q ' 

(4 9) 

,d , d . dp _ 

(4 10) 

,d d . dp _ 

Conversion of this system to cylindrical coordinates (4 8)- (4. 10) takes the 
form (4 2) where 


u = 


V 
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R = 


T = 


/ M cos 0 

COS 0 

sm0 \ 

I cos0 

M cos0 

0 

\ sin0 

0 

M cos 6 J 

— Msin0 

— sm0 

cos 6 

— sm 0 

— M smO 

0 

COS 0 

0 

-M sin0 


Calculation of g(0) for these equations yield. 

1 


(4.11) 


9 ( 0 ) = 


vl — M 2 sin 2 9 + M cos 9 
and the matrix C has the form 


C = 


1 - MQ -Q - R 

-Q 1 -MQ 0 

-iJ 0 1 - MQ 


where 


Q = g cos 9 — g f sin 0, 


R = g sm 0 + g cos 9 
Solutions of (4.6) are given by 


1 


1*MQ 


ao = ho(0) 

V 

1 -AfQ 

and the solutions of (4 7) are given by: 


(4.12) ai = hi (9) 



+ ho(6) 


= ho(0) 




+ h' o (0) 



Here ho(0) and hi(0) are arbitrary functions of 9 and the coefficients b x and 
c t are given by 


&2 = 


(cos# + Mr 2 cos0)/2 + Mr f 2 sin# 
1- MQ 


&3 


(sin 0 + Mrz cos0)/2 + Mr' z sin 9 
1 - MQ 
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sin# + Mr 2 sin 9 
° 2 ~ 1 -MQ 


— cos# 4 Mr 3 sin0 
C3 " 1 - MQ 

Collection of these results in the asymptotic expansions yields (to 0(r 5 / 2 )): 

(4.13) 

ko/o j hifi 

P - r l/2 + r 3/2 ’ 

(4.14) 

. ^0/l , . ^OA 

u — r 2 p 4 b 2 4- ^z {2 C2 ’ 

(4.15) 

1 ^ 0/1 u . * 0/1 
V = T3P+ r3/2 63 + r 3 / 2 C 3 

Differentiating u and v with respect to t and using the result /o = /( 
to 0(r -5 / 2 ), we have 

(4 16) 

p. Kfo 

Ut = T 2 pt 4 b 2 4 3/2 c 2 , 
r r 3 / 2 

(4 17) 

P Kfo 

vt = r 3 pt+ r 3 + c 3 

r 

Finally, noting the term involving can be eliminated from the last two 

equations, we have 

(4.18) 

V 

a p t + P u t + 7 v t = -S, 

T 


where 


a = 7*2 c 3 - r 3 c 2 , 


P = -c 2 , 

7 = C3, 


£ = 63 C 2 — &2 C 3 

The higher order conditions are obtained in a similar manner In fact, one 
can show that the next order condition is of the form 

(ap + P tt + 7 v) tt = —6 + ^re, 

T T* 


(4.19) 
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which is accurate to 0(r“ 7 / 2 ). Note that these conditions do not involve 
any spatial derivatives As such they are ideal for rectangular domains 
where typically one has to pay special attention to the corners, particularly 
when high order numerical schemes are used These conditions correspond 
to the primary acoustic boundary condition. In addition, one must impose 
at inflow boundaries, a vorticity condition For nomsentropic flows, in 
addition to the vorticity, entropy must also be specified at inflow At 
inflow the y momentum equation and zero vorticity condition yield an 
exact relation 

(4 20) v t 4- (p+ Mu) y = 0. 

5. A Model Problem. In this section we begin with the linearized 
Euler equations with mean velocity convection. The scaled form of these 
equations are identical to the one that we used to derive the conditions, 
except they contain forcing terms that characterize a driving source They 
are. 

( 51 ) 

(5 - 2) ^ +M h )%+ % =9i{x ' v ' t) ' 

( 5 - 3 ) 

where gl and g2 model a Gaussian momentum source, which both oscillates 
sinusoidally and decays algebraically in time. Typical examples include, 
a quadrupole sound distribution. Here gi and g 2 sire the gradient of a 
potential <f> Such a function is given by 

<j> = A(t)e~ aR * cos(26) 

where tan# = R = \/(x - xo) 2 + (y — Jto) 2 , (xo,yo) is the location 

of the source, a is a positive constant, and A(t) is the amplitude and a 
function of t. (In the numerical experiments A = sin 27rt/(l + 1 2 ) ) 

In a sample computation which was computed for a time length of 
100 periods (22415 time steps), the solution obtained with the second or- 
der conditions was compared with the exact solution, a solution obtained 
by setting incoming characteristic variables to zero, and one obtained us- 
ing the first order condition. The exact solution was computed in a large 
domain in which, within the time of computations, the waves could not 
reflect off the artificial boundaries and return to the small domain The 
maximum error m pressure calculations observed for the characteristic con- 
ditions was 10.3%, for the first order asymptotic condition it was 3.3%, and 
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for the second order condition the error was 1.3%; indicating the expected 
improvement. In Figures 5.1 and 5.2, the exact solution for the pressure 
is given after 5 periods and 10 periods of time respectively. Subsequent 
pairs of figures (5 3-5.4, 5 5-5 6 and 5 7-5.8) indicate the solution at these 
times for the characteristic boundary condition, the first order asymptotic 
condition, and the second order asymptotic conditions respectively At 10 
periods the errors are visible in the first two cases and their orders of the 
magnitude indeed are as indicated above. Clearly the higher order condi- 
tion improved the results 




FlG. 5.2. Exact Solution at 10 periods of time 
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ion with characteristics based boundary condition t = 10 periods 




FIG. 5.6. Solution with first asymptotic boundary condition t — 10 periods 
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Fig. 5.8. Solution with second asymptotic boundary condition t = 10 periods 
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